The strong and weak formulations for our L-shaped heat transport problem solved by the classical finite element method are the same as for the isogeometric method. \( a(u,v)=l(v) \\ a(u,v) = \int_{\Omega} \left(\frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{ \partial v}{\partial y } \right)v(x,y) dxdy \\ l(v) = \int_{ \partial \Omega_N } g(x,y) v(x,y) dS +\int_{\Omega } f(x,y) v(x,y) dxdy \)
In our example problem of heat flow in an L-shaped area, we assume \( f(x,y)=0 \) (i.e. no heat sources).
On the "outer" part of the shore which we call the Neumann boundary \( \Gamma_N = \{(x,y):x=-1 \cup x=1 \cup y=1 \cup y=-1 \} \)
we assume a known speed of temperature change.
\( \frac{\partial u(x,y)}{\partial n }=g(x,y) \textrm{ dla } (x,y) in \Gamma_N \)
In order to define the Neumann boundary condition, we introduce a polar coordinate system centered at (0,0) and the abscissa lying on the OX axis of the Cartesian coordinate system.
In our example of heat flow in an L-shaped area, we assume the exemplary heat flow velocity function defined in the polar coordinate system
\( R \times (0,2\pi) \ni (r,\theta)\rightarrow g(r,\theta)=r^{2/3 }sin^{2/3 }(\theta+\pi/2) \)
where \( r=(x^2+y^2)^{0.5 } \) is the distance from the origin of the coordinate system, \( \theta = atan2(y,x) \) is the angle between the point \( (x,y) \) and the axis \( x \). So
\( g(r,\theta)=r^{2/3 }sin^{2/3 }(\theta+\pi/2)=(x^2+y^2)^{1/3 }sin^{2/3 }(atan2(y,x)+\pi/2)=g(x,y) \).
On the "inner" part of the shore which we call the Dirichlet shore \( \Gamma_D = \partial \Omega \ \Gamma_N \) we assume a temperature value of zero
\( u(x,y)=0 \textrm{ dla } (x,y) in \Gamma_D \)
An example of a mesh made of triangular elements and the first degree Lagrange polynomials stretched over it is shown in Fig. 1-Fig. 2.

Let us assume that we use the first-degree Lagrange polynomials to approximate the solution. In the case of using linear basis functions, they are related to the vertices of triangular elements (in the sense that their maxima are above the vertices and are defined on the adjacent triangular elements).
If, on the other hand, one were to apply higher degree Lagrange polynomials, they would also be associated with the edges and interiors of the elements. Details are provided in the modules on how to define the basis functions.
In the classical finite element method, our mesh and the resulting matrix of a system of linear equations do not have a regular shape, as is the case with a group of elements in the isogeometric finite element method. Therefore, a global numbering of basis functions related to the vertices of triangular elements should be generated. Each triangular element in turn has its own local vertex numbering. It is necessary to develop a map from global numbering to local numbering of each triangular element. It is drawn in Fig. 3.

Then over the reference triangular element \( \hat{E}:(0,1)-(0,0)-(1,0) \) we introduce three local base functions
\( \hat{\psi}^1(\xi,\eta)=(1-\xi)\eta \\ \hat{\psi}^2(\xi,\eta)=(1-\xi)(1-\eta) \\ \hat{\psi}^3(\xi,\eta)=\xi(1-\eta) \)
These functions are mapped to basis functions spanning over individual elements. For this purpose, we define maps that map the model element to individual elements
\( \hat{E} \ni (\xi,\eta) \rightarrow x_{E1}(\xi,\eta)=(\xi-1,\eta)=(x,y)\in E_1 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E2}(\xi,\eta)=(-\xi,-\eta+1)=(x,y)\in E_2 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E3}(\xi,\eta)=(\xi,\eta)=(x,y)\ni E_3 \\ \hat{E} \in (\xi,\eta) \rightarrow x_{E4}(\xi,\eta)=(\xi+1,\eta+1)=(x,y)\in E_4 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E5}(\xi,\eta)=(\xi,\eta-1)=(x,y)\in E_5 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E6}(\xi,\eta)=(1-\xi,-\eta)=(x,y)\in E_6 \)
Then we calculate the inverse maps
\( E_1 \ni (x,y) \rightarrow x_{E1}^{-1}(x,y)=(x+1,y)=(\xi,\eta)\in \hat{E} \\ E_2 \ni (x,y) \rightarrow x_{E2}^{-1}(x,y)=(-x,-y+1)=(\xi,\eta)\in \hat{E} \\ E_3 \ni (x,y) \rightarrow x_{E3}^{-1}(x,y)=(x,y)=(\xi,\eta)\in \hat{E} \\ E_4 \ni (x,y) \rightarrow x_{E4}^{-1}(x,y)=(x-1,y-1)=(\xi,\eta) \in \hat{E} \\ E_5 \ni (x,y) \rightarrow x_{E5}^{-1}(x,y)=(x,y+1)=(\xi,\eta) \in \hat{E} \\ E_6 \ni (x,y) \rightarrow x_{E6}^{-1}(x,y)=(-x+1,-y)=(\xi,\eta)\in \hat{E} \)
Having mapping formulas between the template element and individual elements, we can define basis functions on the first element
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E1}^{-1}(x,y)\right)=\hat{\psi}^1\left(x+1,y \right)=(1-(x+1))y=-xy \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E1}^{-1}(x,y)\right)=\hat{\psi}^2\left(x+1,y \right)=(1-(x+1))(1-y)=x(y-1) \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E1}^{-1}(x,y)\right)=\hat{\psi}^3\left(x+1,y \right)=(x+1)(1-y) \)
(note that the first local function takes the value 1 at point (-1.1), the second local function takes the value 1 at point (-1.0), and the third local function takes the value 1 at point (0.0)), on the second element
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E2}^{-1}(x,y)\right)=\hat{\psi}^1\left(-x,-y+1 \right)=(1+x)(-y+1)=(1+x)(1-y) \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E2}^{-1}(x,y)\right)=\hat{\psi}^2\left(-x,-y+1 \right)=(1+x)(1-(-y+1))=(1+x)y \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E2}^{-1}(x,y)\right)=\hat{\psi}^3\left(-x,-y+1 \right)=-x(1-(-y+1))=-xy \)
(note that the first local function is 1 at (0,0), the second local function is 1 at (0,1), and the third local function is 1 at (-1,1)), on the third element
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E3}^{-1}(x,y)\right)=\hat{\psi}^1\left(x,y \right)=(1-x)y \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E3}^{-1}(x,y)\right)=\hat{\psi}^2\left(x,y \right)=(1-x)(1-y) \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E3}^{-1}(x,y)\right)=\hat{\psi}^3\left(x,y \right)=x(1-y) \)
(note that the first local function is 1 at (0,1), the second local function is 1 at (0,0), and the third local function is 1 at (1,0)), on the fourth element
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E4}^{-1}(x,y)\right)=\hat{\psi}^1\left(x-1,y-1 \right)=(1-(x-1))(y-1)=x(1-y) \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E4}^{-1}(x,y)\right)=\hat{\psi}^2\left(x-1,y-1 \right)=(1-(x-1))(1-(y-1))=xy \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E4}^{-1}(x,y)\right)=\hat{\psi}^3\left(x-1,y-1 \right)=(x-1)(1-(y-1))=(1-x)y \)
(note that the first local function takes the value 1 at point (1.0), the second local function takes the value 1 at point (1.1), and the third local function takes the value 1 at point (0.1)), on the fifth element
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E5}^{-1}(x,y)\right)=\hat{\psi}^1\left(x,y+1 \right)=(1-x)(y+1) \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E5}^{-1}(x,y)\right)=\hat{\psi}^2\left(x,y+1 \right)=(1-x)(1-(y+1))=(x-1)y \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E5}^{-1}(x,y)\right)=\hat{\psi}^3\left(x,y+1)\right)=x(1-(y+1))=-xy \)
(note that the first local function is 1 at (0,0), the second local function is 1 at (0, -1), and the third local function is 1 at (1, -1)),
and on the sixth element
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E6}^{-1}(x,y)\right)=\hat{\psi}^1\left(-x+1,-y \right)=(1-(-x+1))(-y)=-xy \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E6}^{-1}(x,y)\right)=\hat{\psi}^2\left(-x+1,-y \right)=(1-(-x+1))(1-(-y))=x(1-y) \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E6}^{-1}(x,y)\right)=\hat{\psi}^3\left(-x+1,-y)\right)=(1-x)(1+y) \)
(Note that the first local function is 1 at (1, -1), the second local function is 1 at (1,0), and the third local function is 1 at (0,0)).
Functions on local elements are aggregated into global base functions \( e^1,e^2,e^3,e^6,e^8 \) according to the mapping shown in Fig. 3.
In other words
\( e_1=\psi^3_{E1}+\psi^2_{E2} \\ e_2=\psi^2_{E2}+\psi^3_{E3}+\psi^2_{E4} \\ e_3=\psi^1_{E4}\\ e_6 =\psi^2_{E3}+\psi^3_{E4}+\psi^1_{E6} \\ e_8=\psi^2_{E5}+\psi^3_{E6} \)
We approximate our solution (temperature distribution) by a combination of aggregated Lagrange polynomials
\( u(x,y)=\sum_{i=1,...,8 } u_i e_i(x,y) \).
Because of the use of the zero Dirichlet condition on the "inner" edge, we can immediately assume that the Lagrange polynomial coefficients related to the inner edge are equal to zero, i.e.
\( u_4=u_5=u_7=0 \), so we didn't even derive the formulas for base functions \( e_4(x,y),e_5(x,y),e_7(x,y) \), because they are reset by a zero Dirichlet condition (forced by zeroing the coefficients).
Generation of irregular global numbering of basis functions and irregular mapping of local element functions to global functions is a typical feature of the application of the finite element method on irregular meshes, for example with the use of Lagrange polynomials. In general, we also need to generate the global numbering and local numbering of the vertices of the triangular elements, and generate the formulas for the basis functions on all triangular elements